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Abstract 

The anisotropy of wood within the radial-tangential (RT) growth plane has a major in- 
fluence on the cracking behavior perpendicular to grain. Within the scope of this work, a 
two-dimensional discrete element model is developed, consisting of beam elements for the 
representation of the micro structure of wood. Molecular dynamics simulation is used to 
follow the time evolution of the model system during the damage evolution in the RT plane 
under various loading conditions. It is shown that the results are in good agreement with 
experiments on spruce wood, and that the presented discrete element approach is applica- 
ble for detailed studies of the dependence of the micro structure on mesoscopic damage 
mechanism and dynamics of crack propagation in micro structured and cellular materials 
like wood. 
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1 Introduction 



Wood is probably the most ancient structural material in the world, widely used 
today in many species, for all kinds of purposes, leading to a world production of 
roughly 10 9 tonnes, which is about even to the one of iron and steel [1]. It's ap- 
plication in modern civil constructions of large dimensions like stadium roofs or 
long- span bridges calls for a good understanding of design properties like mod- 
uli, crushing strength and toughness. Wood is a natural fiber composite, which is 
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strongly orthotropic due to its internal material structure. The elemental material is 
a cellulose cell wall material with equal properties for different species [1]. Never- 
theless, macroscopic properties vary highly within the tree, among different species 
and along loading directions. Explanations for these extreme differences are found 
on the microstructural level in the cellular structure. Several internal length scales 
can be found, that are relevant for the damage evolution, including many defects. 
Depending on the observation scale and position, wood shows different physical 
and geometrical properties. The size scales are typically subdivided into hierar- 
chical levels like the atomic, micro, meso and macro scale. The material can be 
regarded as a cellular structure with hexagonal shaped wood fibers arranged in par- 
allel on the micro scale (comp. Fig.l). The dominating structure on meso scale in 
the RT plane, perpendicular to the tree axis, are annual rings with an alternating 
early- and latewood that differ in fiber shape and thickness of the cell walls. For 
the understanding of damage in wood, damage mechanisms on the micro and meso 
scale are of special interest. 




Figure 1. Relevant scales in wood (Micro scale picture from [8]) 

Discrete Element Methods (DEM) are a fast emerging computational method de- 
signed to solve problems with gross discontinuous materials and geometrical be- 
havior. Due to the creation and continuous motion of evolving crack surfaces, frac- 
ture of material is difficult to handle numerically. Continuum models have problems 
to account for the discrete nature of material failure. One approach is a smeared 
crack approach, e.g. with a micro-plane damage model [2], where the influence of 
the damage is considered in terms of the material properties. The other approach 
is to continuously change the model topology, as cracks evolve [3]. Alternatively, 
discrete models like lattice dynamics models can be used to simulate fracture [4]. 
A common type of question around DEM is related to the possibility of assign- 
ing continuum properties to a given lattice structure, especially avoiding artifacts 
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e.q. due to a regular micro structure of a model. Contrary, for the representation 
of wood, a regular microstructure like a hexagonal lattice is especially suited for 
naturally taking the anisotropy of wood in the RT plane into account. 



In the present paper fracture of wood was studied with a combined beam-particle 
lattice model on the micro- and mesoscopic scale. The mesoscopic structure of 
annual rings is projected on the microscopic model and generic failure mechanisms, 
observed in the failure of spruce wood in the RT plane, are implemented into the 
model. The paper is organized as follows: Section 2 gives a detailed description 
on experimental in-situ studies of the fracture processes in clear softwood under 
various loading conditions perpendicular to grain. Experiments are focused on the 
micro-meso- structural scale, using a confocal laser scanning microscope (CLSM). 
Identified generic damage mechanisms are the rupture of earlywood cell walls for 
crack propagation in tangential direction and debonding of the interface of wood 
fibers for crack propagation in radial direction. Section 3 provides an outline of 
the theoretical background of the discrete element model, for the micro-structure 
of wood. After the model description, test simulations and damage simulations are 
presented and compared to the experimental findings in Section 4. 



2 Damage in wood in the RT-plane 



Damage basically initiates on the atomic scale and reaches relevance for larger 
scales like the micro or meso scale while it propagates, leading to global failure 
when reaching the macro scale. For the prediction of the damage evolution in ma- 
terials scale dependent knowledge on their relevance for the damage evolution is 
indispensable and intimately connected to its characteristic structure. The typical 
observation scales of wood are shown in Fig.l. Looking at the sample of wood suf- 
ficiently distant from the center of the tree, three orthogonal planes of symmetry 
are found: radial (R), tangential (T) and longitudinal (L). Stiffness and strength 
are greatest in longitudinal direction. In real constructions load perpendicular to 
grain is often present and damage consequently. Macroscopic material values for 
spruce wood are given in the table [3,5]. Even though wood, as we use it, is just the 
dead corpus of a living tree, it is the result of its process of formation. Therefore, 
deviations within one board of wood can be in the range of 25% [5], since during 
the growth process, trees have to react to environmental conditions by structural 
adaptation on all scales. 
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Typical properties for spruce wood [3]: 

Young's / shear moduli [MPa] I Poisson's ratios: 
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On the meso scale, wood is a cellular solid with mainly cells of hexagonal cross- 
section (also called tracheids in the case of soft-wood). The characteristic annual 
rings in soft-wood consist of alternating circumferential bands of thick- and thin- 
walled tracheids. Fig. 2 shows new measurements of the relative density distribu- 
tion p*/p c , based on CLSM pictures for an annual ring with the cell wall density 
p c . Other features of the meso scale are rays and sap channels, that are important 
for the water balance and growth of trees but are of minor importance for the ma- 
terial representation of soft-wood [1,3,6]. Of course, structural differences between 
species are significant. Micro mechanical finite element models reveal the cell den- 
sity to be the governing parameter of the macroscopic mechanical properties of 
wood within the RT plane [6]. 




Figure 2. Optically measured relative density, cell aspect ratio and microstructure in one 
annual ring 

Fibers are glued together by lignin, leading to a natural fiber composite of less 
than 5% matrix volume fraction. On the finer micro scale, wood is again a fiber- 
reinforced composite, since the cell walls are made up of fibres of cellulose, which 
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is a crystalline polymer with sort amorphous regions, embedded in an amorphous 
hemi-cellulose and lignin matrix [6]. These micro-fibrils are oriented in layers of 
different orientation, depending on the time of formation [1,7,8]. 




Figure 3. Specimen with dimensions in mm and different investigated configurations of 
load to growth directions 



Cell wall properties for softwood [1]: 



cell wall density p c 

Young's modulus(T) E T 

Shear modulus G T 

Yield strength (T) it c 
Volume fraction 



[kg/m 3 ] 

[GPa] 

[GPa] 

[MPa] 

[%} 



1500 
10 

2,6 
50 

85 - 95 



Tracheid size longitudinal [mm] 

tangential size It 

radial size l R [fim] 

wall thickness t [fim] 



2,5 
30 

12 - 40 

1,6-5 



7,0 



Even though different kinds of soft wood have a large difference in relative density, 
the cell wall properties for all species are pretty much identical. This is reasonable, 
since in nature, optimization of structural materials is only possible in the frame- 
work of feasible biological processes, which are identical for all trees. Therefore 
optimization of microstructure and shape is indispensable for survival. To under- 
stand the mechanical and breakdown behavior of wood, knowledge on the material- 
structure composition of the high strength elements of wood - the cells is necessary. 



2.1 Microscopic in-situ experiments 



The damage evolution in wood on the micro-meso structural level in the RT-plane 
was studied under tension. All test specimen are clear spruce wood, cut from one 
board with a mean density of 500kg /m 3 and a moisture content of 12%. The annual 
ring width dj r varied between 1mm and 2.5mm for different specimen, due to 
different locations within the board cross section. Details on specimen preparation 
and experimental setup are published in [9,10]. 

In-situ tests are made using a Confocal Laser-Scanning Microscope (CLSM), there- 



5 



"b.<> 0,1 <12 03 0.4 0.6 O.J 0,8 

TR RT RT45° displacement [mm] 

Figure 4. Comparison of crack patterns for different configurations (top) and load displace- 
ment curves for TR, RT and i?T45°-configuration (bottom). 

fore intensities on the pictures represent height values from a focus plane. 



2.2 Micro- and mesoscopic degradation of wood 



Different orientations lead to completely different damage evolutions. If wood is 
loaded in tangential direction (Ti?-specimen (Fig. 4a)) the crack propagates straight 
in radial direction with minor deviations and smooth crack surfaces. For the RT 
configuration (comp. Fig. 46) the crack propagates within the earlywood parallel 
to the annual ring and perpendicular to the loading direction. The crack surface is 
rough with edges and steps and the crack propagates mostly in an instable man- 
ner. Compared to pure RT an TR crack configurations a zig-zag behavior is ob- 
served for orientations under 45° (i?T45° Fig.4c). Initiated at the notch, the crack 
first propagates in radial direction with smooth crack surfaces. While the crack ap- 
proaches the early /latewood transition, it is deflected from the initial direction until 
it propagates along the transition edge. The final failure occurs, when half of the 
specimen is damaged, by a sudden crack formation in radial direction. 

Microscopic observations reveal different damage mechanisms. Tension stress in 
radial direction is mainly released by rupturing the cell walls, leading to rough and 
stepped surfaces (comp. FigAbottom insert 6) while for tension stress in tangential 
direction cells are generally intact, but inter-cellular damage (fiber debonding) is 
the dominant damage mechanism (comp. FigAbottom insert a). For the RTA5° 
configuration we find a combination of inter- and intracellular damage (comp. 
FigAbottom insert c), depending on the direction of the crack propagation. Fur- 
ther details on the damage evolution under tension can be found in [9,10]. 
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3 A discrete element model for wood 



In order to study the damage evolution of wood, a two-dimensional network model 
is worked out. Molecular Dynamics (MD) simulation is used to follow the motion 
of each cell by solving Newton's equations of motion. In the present study we use 
a A th order Gear Predictor Corrector scheme [11]. A general overview of MD 
simulations applied to composite materials can be found in [4]. 




Figure 5. Degrees of freedom and internal loads on a single beam element and hexagonal 
lattices. 

With this method already small cracks are sharply defined, with the possibility to 
simulate simultaneously a conglomerate of cracks within rather large lattice sizes 
[4]. The fundamental advantage of the lattice model used in this investigation is 
due to its simplicity, giving direct access and possible physical interpretation to 
each step of the algorithm. Consequently, one can modify the rules for features 
of interest, like the characteristic properties of size, strength or force in a rather 
straightforward and transparent way. 

The model is composed in four major steps, namely, a) the implementation of 
the microscopic structure, b) the determination of the constitutive behavior, c) the 
breaking of the model and finally d) the construction of the test specimen. 

(a) Micro structure: The basic structural element used is a two-dimensional TlM- 
OSHENKO combined spring-beam element for small displacements, in other words 
elements, that can elongate like springs and bend like beams at the same time. An 
overview over network models is given in [12]. The nodes with mass % — 1...N 
of a two-dimensional lattice with N nodes are connected via these elastic elements. 
An element connecting the positions i and j, has the cross section A l i and the geo- 
metrical moment of inertia J' 3 ' of the beam for flexion. The length l l i of all beams 
is equal for this case. The elastic behavior is also governed by the two material de- 
pendent constants for the Young's- and shear moduli of the beam E b and Gb- For 
a beam between sites % and j the normal, shear and bending flexibilities o^, b^, 
are given by: 



Each node % has in the global coordinate system the position x^ yi and the rotation 
9i and all the forces and moments of the beams, connected to i act on it, allowing 
only motion in the observation plane (comp. Fig. 5a). In the local coordinate system 
of the beam, three continuous degrees of freedom are assigned to both ending points 
i, j of the beam, which are for site % the displacements Ui and Vi and the bending 
angle The longitudinal force acting at site % is 

F b,x = a ij( u i - u i)i witn «<j = IMjj ( 2 ) 

the shear force is 

Kv = Pi&j-Vi)-^(<l>i + <l>i), with = 1/(^ + ^-/12), (3) 
and the flexural torque at site % is 

K,z = "7^0 '■■l, l O J )-<); J lj l (0 J - O,). With % = /^(fetf/Ctf+1/3), (4) 

consequently resembling the Timoshenko beam theory. Quantity 6^ was chosen 
to be b tj = 2a lj , resulting in a PoiSSON's ratio of v b = 0. The model is derived from 
a model for geomaterials, presented in [13,14]. For simplification purposes, the 
beams behave in a time independent way, but time dependent beam properties are 
easy to introduce. In our simulation, the nodes are the smallest particles interacting 
elastically with each other. 

The construction of the final model is done in several steps shown in Fig. 7. First 
a hexagonal beam lattice is constructed. As second step, the meso structure of the 
annual rings is projected on the fundamental lattice with the possibility to vary the 
distance d 3r and angle a of the annual rings with respect to the specimen axis, which 
is identical to the loading direction. The different material properties of an annual 
ring are represented by alternating the effective density p(x' lj ") = p* according to 
the function 

P* = Pmin (l + 2 • axf 6 • e«?") (5) 

representing an annual ring. If the sample of wood is cut distant from the center of 
the tree, the curvature of annual rings can be neglected, leading to orthotropic mate- 
rial properties. Eq.5 is used in order to obtain a steady density function. Dimension 
x l J stands for the distance of the middle point of the beam from the annual ring 
interface, connecting sites i and j in relative coordinates of the R — T— material 
coordinate system (compare Fig. 7). The parameters a, 6, c, are selected in a way, 
that the properties scale according to the density between an annual ring (comp. 
Fig. 2). The thickness for the beam connecting site % with j depends on the lattice 
type and parameters chosen. 

The required relationship between the cell wall thickness Wij, the edge lengths h 
and /, and the relative density p* / p c for hexagonal lattices with corrections for 
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Figure 6. Construction of the model with ground lattice, meso structure projection and 
cutting of the specimen 



overlapping areas is 



w ij (p*) = A- \(h + 2l)- 




+ 2Z) 2 -cosfl(2/i + Zsinfl)-p- 

Ap t 




(6) 



with A = 1 + sin 6 + cos 0. 

(b) Constitutive behavior: Since the model is of two-dimensional nature, we an- 
alyze its response to load applied in the R — T— plane. Even though all elements 
are considered to be isotropic, anisotropy arises due to their arrangement and the 
macroscopic moduli E* R) E^, G* RT , v RT (since E R u^ R = E^v^ R is assumed) and 
plateau values for the compressive stresses a* R and er£ if compressive failure is 
dominant, are required to describe the mechanical behavior of the composite. The 
derivation for the relations assuming linear elastic deformations can be found in 



For regular hexagonal lattices with the parameters h, I, 6 and the beam thickness w 
the Young's moduli and E* R and the shear modulus G* RT axe calculated with 



[1]. 
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the properties of the cell wall E c as 



E* R _ (w\ 3 cos9 E?p _ fw\ 3 (h/l + sini 



E c \IJ (h/l + sin 9) sin 2 6' E c \lj cos 3 9 ' 
g *rt _ f w \ 3 {h/l + ame) „ cos 2 # () 



E c \lj (h/l) 2 (l + 2h/l)cos9' Urt (h/l + sm9)sm9' 

With these expressions it is possible to estimate the correct values for the mi- 
crostructure iteratively, since a change of any property also effects the other prop- 
erties through the density function Eq.5. Good micro structural values for spruce 
wood were found to be E c = llGPa, h = 1.5/ and 9 = 12°, but of course these 
values depend on the number of wood fibers, that form a annual ring, leading to 
cell wall thicknesses 2.6/im < Wij < 10.3/xm. These values do not fit the observed 
microstructure, since Eqs.7 are only valid for regular hexagonal lattices. Solutions 
for the real model are presented in Section 4. 



(c) Breaking of the material: Stresses inside the system can be released by the frac- 
ture of the beam elements, either under tension or under bending. In the framework 
of Discrete Element Methods, the complicated crack-crack interaction is naturally 
taken into account. If the deformation of an element due to the total forces acting 
on it exceeds its breaking thresholds (damage thresholds comp. Fig. 6), its stiffness 
is abruptly reduced to zero, resulting in a load redistribution to neighboring beams 
in the next iteration steps. If after some iterations due to the load redistribution to 
the neighboring beams they exceed their threshold value due to the additional load, 
they fail too, giving rise to crack growth. Alternatively, so called continuous dam- 
age laws with incremental stiffness reduction, once the actual damage threshold is 
exceeded, can be used [15]. Since we have no particular basis for a topological dis- 
order or cell property disorder as well, we make a computational simplification by 
assuming the disorder in beam properties to adequately account for all the relevant 
disorder present in the represented material region. Disorder is introduced to the 
model by a two parameter Weibull distribution of breaking thresholds e d at the be- 
ginning of the simulation and is kept constant during the fracture process (quenched 
disorder) as 

P(eo) = l-exp((-^) m ), (8) 

with the Weibull modulus or shape parameter m = 6 and the scale parameter e = 
0.0055. 



The two experimentally identified stress release mechanisms, cell wall breaking and 
peeling are represented in different ways. The peeling mechanism (comp. Fig. 6) is 
activated, when the force on a node in tangential direction exceeds a force threshold 
value. This value is calculated by multiplying a critical debonding stress pi and 
the area A affected by the debonding, since the matrix material is the same for 
the whole wood material (comp. Fig. 6). Contrary, the cell wall breaking depends 
on the cell wall thickness. The breaking criteria always use the same disorder for 
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Figure 7. Representation of different damage mechanisms 

all damage mechanisms for simplicity. p 2 is a parameter to take the quantitative 
relation of the critical values for tension and bending into account. 



(d) Construction of the specimen: Finally, the simulation model is extracted out 
of the basic hexagonal lattice under the angle [3 by deleting those elements, which 
do not fit in the specimen geometry given by Fig. 3. As a last step, the border ele- 
ments are identified and connected to the two loading points via elastic beams. This 
way, the specimen is free to rotate as observed in the experiments, when the crack 
propagates. 



4 Computer simulations 



After the model construction, computer simulations for the constitutive behavior 
of the material are performed. This is done by calculating the time evolution of 
the particle assembly by solving the equation of motion of the nodes in terms of 
molecular dynamics, 

rm% = Pu i = l,...,N (9) 
where N denotes the number of nodes, is the mass and Fj is the total force acting 
on node %. A A th order Predictor Corrector algorithm is used in the simulations to 
solve numerically the second order differential equation system (Eq.9) [11]. After 
each integration step the breaking condition is evaluated for all the intact beams. 
Those beams which fulfill this condition are removed from the further calculations. 
The discrete element method gives the possibility to monitor the development of 
the microscopic damage in the specimen. 

Material properties are measured numerically on the macroscopic scale by com- 
puter simulations on the bulk model. All beams are considered to be linear elastic up 
to failure, so no plastic or time dependent creep behaviour is considered in this stage 
of the analysis, but basically the elements can easily be replaced by other kinds of 
rheological elements. Fig. 8 shows the constitutive behavior of a model block with 
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Figure 8. Stress-strain curve for a sample of lxlmm with annual ring distance of 1mm 

different material orientations. Values are transformed into non-dimensional form 
for comparison, using characteristic values of single beams . The increase of the 
radial Young's modulus E R is due to the alignment of the beams in loading direc- 
tion. 

The Young's modulus is calculated via the energy e pot stored in the beams as 
E R / T = 2e pot /V ■ e~ 2 for the volume V, and is in accordance with values for 
spruce wood. 

As a next step, breaking is switched on and the three configurations RT/TR/A5°RT 
are calculated. Activated crack mechanisms and crack evolution are in agreement 
with the experimental observations. 



Figure 9. Microstructure of damage for the three different configurations RT; TR; 45° RT 

As a further application, we simulate the permeation of a circular object (e.g. a nail) 
in the center of our model by defining a repulsive contact force. During the simu- 
lation we constantly increase the diameter of the nail, leading to damage. Contact 
forces are only calculated for the interaction of the nail with the cell structure and 
not within the cell structure itself. We observe crack growth mainly in the tangen- 
tial direction. A detailed study on the compressive failure of cellular materials with 
the discrete element method will be published in a follow-up paper. 
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Figure 10. Snapshot of the damage for the permeation of a nail. 



5 Discussion 



A disordered hexagonal beam network model has been introduced to study the dam- 
age development in the natural fiber composite wood. The main advantage of our 
modeling is, that it naturally accounts for the complicated local stress fields formed 
around failed regions in wood. Furthermore, it captures the gradual activation of 
the two relevant damage mechanisms and their interaction during the fracture pro- 
cess. We have demonstrated that our discrete element model is capable to provide a 
deep insight into the damage process occurring perpendicular to grain under grad- 
ual loading of wood. Experimental results on the micro structure of damage and on 
the loss of stiffness are in good agreement with numerical simulations. With our 
model even dynamics of fracture during the processing of wood can be simulated, 
as demonstrated on the damage development caused by the permeation of a nail. 



An extension of the model to the description of time-dependent failure processes is 
straight forward, since time-dependent rheological elements can easily be inserted 
on the microstructural level. The capabilities of our model are not limited to wood, 
it can be applied to study dynamic damage processes of other cellular or micro- 
structured materials under varying loading cases including also thermal or chemical 
degradation processes resulting in complex constitutive governing equations of the 
single elements. Studies in this direction are in progress. 
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